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Abstract. The Continnum Dislocation Dynamics (CDD) theory and the Discrete 
Dislocation Dynamics (DDD) method are compared based on concise mathematical 
formulations of the coarse graining of discrete data. A numerical tool for converting 
from a discrete to a continuum representation of a given dislocation conhguration is 
developed, which allows to directly compare both simulation approaches based on 
continuum quantities (e.g. scalar density, geometrically necessary densities, mean 
curvature). Investigating the evolution of selected dislocation configurations within 
analytically given velocity fields for both DDD and CDD reveals that CDD contains a 
surprising number of important microstructural details. 
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1. Introduction 

Predicting and understanding the collective evolution of dislocation ensembles together 
with the resulting mechanical properties of crystalline materials is a long-standing 
goal of microstructure-based plasticity. From a computational perspective, reducing 
dislocation mechanics to a system of partial differential equations amenable to solution 
by standard continuum computational methods would be highly desirable. However, 
the nature of those curved, line-like defects makes the development of a continuum 
theory of dislocations an extremely difficult challenge (^, which today is still not fully 
mastered. One way of enriching macroscopic continuum descriptions with additional 
details of smaller length scales (e.g. with details of dislocations) are concurrent as 
well as hierarchical multi-scale approaches for bridging the atomistic to the continuum 
domains [2]-|^. These methods aim at establishing a direct link between atomistic details 
of dislocations and macroscopic mechanical properties of materials. The level of detail 
in these approaches is very high, and length scales can be reached that are not accessible 
to one-scale methods as e.g. atomistic simulations. Multi-scale approaches, however, are 
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still computationally very expensive and the information transfer between methods is 
numerically challenging. An alternative description of dislocation ensembles is utilized 
in mesoscopic methods such as the discrete dislocation dynamics method (DDD) [S 15 


In DDD the atomic scale is not resolved, and instead dislocations are represented as 


discrete, linear defects interacting according to the elastic theory of dislocations 16 20 


Because each dislocation is resolved individually, DDD is able to accurately describe 
the motion and interaction of dislocations in a very detailed manner. However, the 
computational cost of DDD scales with the number of dislocations considered and is 
therefore computationally expensive when it comes to large numbers of interacting 
dislocations {large referring e.g. to a density of >10^^m“^ in a volume of >(10pm)^). 
To overcome these limitations, one might seek continuous density-based descriptions as 
alternative to the representation of dislocations as discrete objects. Reducing dislocation 
microstructure to distinct “types” of e.g. geometrically necessary and statistically stored 


dislocations (GNDs and SSDs) allows strain gradient based continuum methods 21 23 


to become independent of the number of interacting dislocations. These approaches 
result in a substantial gain in terms of tractable length and time scales which makes 
them interesting for engineering applications. The main drawback is that only few, 
strongly averaged aspects of systems of dislocations are preserved, and e.g. fluxes 
of dislocations cannot be accounted for at all. Recently emerging continuum models 
of dislocation dynamics (CDD) bridge the gap between these last two approaches by 
combining the efficiency of a continuum theory with the physically sound basis of a 
dislocation based description: ’Kroner-Nye tensor’-related models 24-26 are based 


on fluxes of GNDs, and continuum screw-edge representations are a relatively coarse 
but efficient way of representing dislocation loops 27-29 , while models that consider 


densities of curved lines contain considerably more information 30 35 and recently 


even have been coupled to gradient plasticity models [^. Among others, applications 
with high densities and/or high accumulated plastic strains particularly benefit from 


continuum dislocation dynamics models (e.g. 37 ^). 

Because in a GDD-type continuum framework details of individual dislocations 
are accounted for in an average sense, one may question the physical validity and 
the predictive capability of such an approach. In principle it should be possible to 
benchmark GDD theories through comparison with averaged data extracted from DDD 
simulations. Up to date, contributions into this direction have been presented by A. 
El-Azab and co-workers 42,43 , who investigated ensemble averages and statistical 
properties of dislocation microstructure e.g. in terms of dislocation density and 
orientation distributions. The objective of the present work is to compare GDD 
and DDD directly by analyzing dislocation microstructures that develop in time from 
equivalent initial conditions. Although such a detailed comparison is essential to validate 
assumptions used to derive the particular GDD formulation, no previous attempts into 
this direction exist (to the best of the authors’ knowledge). One of the reasons for 
the lack of such direct comparisons may reside in the formulation of a GDD theory 
itself. In order to be comparable to DDD simulations in an average sense a GDD 
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framework must fulfill a number of criteria: (i) its density measures must be derived by 
well defined statistical averaging steps from systems of discrete dislocations (i.e. in a 
bottom-up approach); (ii) its evolution equations should map an initial set of density 
measures onto a set of temporally evolved densities; (iii) it must be kinematically and 
(iv) dynamically consistent. 

• Kinematic consistency describes the ability to evolve dislocations as curved and 
connected lines with the concomitant line length change during motion in a given 
velocity held. 

• Dynamic consistency, on the other hand, describes the ability to predict the actual 
evolution of a dislocation system under externally applied mechanical load as well 
as mutual interactions between dislocations. 


It is generally not recognized in the literature that kinematic consistency is a necessary 
prerequisite for dynamic consistency and considers only geometrical aspects of averaged 
systems of dislocations. While the problem of dynamic consistency is to date still far 
from being solved, the recently developed (so-called ’simplihed’ or ’integrated’) CDD 
theory by T. Hochrainer and co-workers was rigorously derived to solve 

the problem of kinematic consistency and has already demonstrated its applicability in a 
number of problems 34,35 , including dislocation patterning 41 . The goal of this work 


is to establish a methodology for a direct and detailed kinematic comparison between 
DDD and CDD models and to demonstrate the accuracy of CDD through a number of 
benchmark problems. 


The outline of this paper is as follows: in section we give a brief overview of 
the general concept for validating continuum dislocation microstructure with discrete 
data used in this paper. In section we introduce the mathematical foundations for 
the geometrical representation of curved dislocations lines and for obtaining averaged 
continuum quantities in an analytical as well as in a numerical way. We then briehy 
introduce the evolution equations for CDD and how CDD data can be numerically 
obtained from DDD in section Section investigates numerically the quality of the 
continuum model during time evolution of DDD and CDD microstructures for three 
different benchmark systems. Finally, in section we demonstrate how the DDD model 
can be used for “data mining” in order to estimate the quality of one of the assumptions 
on which the CDD version used in this paper is based. 


2. Outline of the CDD-DDD Validation Approach 

The fundamental idea for validating the microstructure evolution within CDD is 
sketched in Fig. Consider an ensemble of dislocation loops with the same Burgers 
vector b, occupying a domain of size y. Ly x Az (Fig. IT). where the slip system 
normal points in z direction. Each dislocation is located on a different slip plane and 
moves - driven by an external stress held - by glide only; no annihilation, cross slip 
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Figure 1: Sketch of the system (A), coarse graining steps and time evolntion. The 
continnous field (C) is obtained from the DDD conhguration (B). Both systems then 
evolve independently. DDD and CDD result in (D) and (E), respectively, and can then 
directly be compared by converting (D) into continuous helds (F). 


etc are considered. Averaging over the height of a snb-volnme Az we obtain Fig. 
which serves as initial dislocation microstrnctnre for onr investigations. In the following 
sections a systematic method for extracting average helds (e.g. dislocation densities, 
Nye tensor, curvatnre, etc) from systems of discrete loops is developed. This discrete- 
to-continuum conversion method {D2C) is applied to the DDD initial conhguration 
(Fig. 133), from which we obtain average helds (Fig. [^). Those serve as initial values 
for the CDD model. CDD initial structure (Fig. [^) is evolved by time integration 
of the CDD evolntion eqnations, while the initial conhguration of discrete loops is 
independently evolved in time by the DDD eqnation of motion. This resnlts in the 
DDD and CDD microstrnctnres as shown in Fig. EP and E. At any snbseqnent point 
of time, relevant state qnantities of CDD, Fig. [pH, can be compared to corresponding 
DDD quantities as obtained from the D2C conversion. Fig. EP- 

3. Conversion of discrete dislocation lines into continuons fields (D2C) 

In what follows the discrete-to-continnnm {D2C) approach is described which is nsed 
to systematically average the geometric properties of discrete dislocation lines to obtain 
CDD state qnantities. To diherentiate between variables for discrete and continuous 
measures the superscript ’d’ is used for discrete quantities. The snperscript is dropped 
for a continuous quantity or when the context is clear. This section starts by introducing 
the geometrical representation of discrete lines and their approximation throngh splines, 
followed by a formal density dehnition. Subsequently, mathematical averaging operators 
are introduced and concisely derived in a discretized form which is snitable for numerical 
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3.1. Discrete dislocation lines and their geometrical approximation 


From a geometrical viewpoint, a dislocation line is an oriented curve c(£) which can be 
parameterized by its arc-length i. The local unit tangent vectors are the first derivative 
of c 


= ( 1 ) 

while the curvature vector is given by the change of orientation per unit arc-length, 




d"c 


di dC 


( 2 ) 


The scalar curvature is the reciprocal of the local radius of curvature R of the dislocation 
line: 


= II 


1 

R’ 


(3) 


Although parameterization by arc-length is the most natural choice for defining 
geometrical properties of a curve, in general the parameter i is only available after 
the curve has been constructed. Thus, for practical purposes, it is often necessary to 
prescribe a curve through a parameter u, say u G [0,1]. In this case unit tangent and 
curvature become 


CM = 


and k'^{u) = 


1 d"c 


J du ” -- V / J2 

where the scalar J and its derivative are given by 


1 dJ dc 
du du 


(4) 


Jiu) = 


d^ 

du 


dc 

du 


and 


dJ 1 dc d^c 
du J du du^ 


(5) 


In this work dislocation lines are discretized into segments, each of which is represented 
by a cubic Hermite spline of the form 


’ 1 000 ' 


Po 

0 

0 

0 


To 

—3 —2'y 3 —7 


Pi 

2 7 —2 7 


Ti 


( 6 ) 


where {Po,To) are position and parametric tangent vectors of the first knot of the 
segment (u = 0 ), (Pi,Ti) are the corresponding quantities for the second knot of the 
segment (u = 1 ), 7 = ||Pi — Po||“, and a is a tension parametei[|j Cubic splines are 
used instead of linear segments in order to have non-vanishing curvature along each 
segment. This type of representation is adopted in the so-called Parametric Dislocation 
Dynamics method (PDD) [^. The numerical studies below are based on the existing 
implementation 


f Parameterization corresponding toQ! = 0,a = 0.5, and a = 1 are called uniform, centripetal, and 
chordal, respectively. 
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As an auxiliary functional for averaging geometrical properties of discrete dislocations, 
we introduce the Dirac delta-operator on a dislocation line c as 

Sa{[»],r) = j 6{c{i)-r)»di, (7) 

0 

where the symbol is a placeholder for an arbitrary line property described by a 
scalar, vectorial or tensorial function of c. Using Q, for each discrete line c we can 
dehne the scalar dislocation density pc, the line direction density (or vector of signed 
GND densities) Qc, the Kroner-Nye tensor olc [^, and the curvature density vector 

af] 

Pc{r) = {.[l],r) = j (5(c(£) - r)di (8) 

0 


Qc{r) = (5c 



La 

jHcW-r)'^de 

0 


(9) 


o^c{r) = 6c 


dc 

— ® h 
di 


,r ] = f 6{c{i) -r)^®bdi 


( 10 ) 


Qcir) = 6c 


^d^c' 


,r]= f 6{c{i)-r)^di. 


( 11 ) 


The later on used scalar, signed curvature density Qc can be obtained by projecting the 
curvature vector on the outwards pointing normal 

La 


qc{r)= [6{c{i)-r)^- l^xn] di, 


dC \ di 


( 12 ) 


where n is the slip plane normal. These discrete measures are well suited for spatial (or 
ensemble) averaging over a number of lines, and we introduce an averaging operator for 
the volume Vj. of size V as 

(•)v',r := (1/V") [ • dV, 


(13) 




where the averaging volume is centered around r. Application of (13) to ([^-(12) leads 
to the dehnition of the corresponding average helds: 

1 


V,r c ' 


di 


'<V,r 


(14) 


§ In the following we assume only one slip system, i.e. all dislocations have the same Burgers vector 
b. From this case we can then construct the multi-slip system situation by appropriate superposition 
(e.g. summation of densities). 
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(15) 



(16) 



(17) 


Therein, C c denotes a section of line c contained inside the averaging volume 


V which is centered at r. Note, that upon volume integration of the total density (14) 


the total dislocation line length contained within the volume is obtained. Analogously, 


integration of the signed curvature density (17) yields the number of closed loops as 
multiple of 27i. The vector of signed GND densities g contains densities of positive 
and negative edge and screw dislocations as later used in the CDD evolution equation 
( [2^ . We remark that the (scalar) geometrically necessary dislocation density pc is the 
norm of g, i.e. pc = g = ||^||, and that the average line direction of the geometrically 
necessary dislocations is given by the unit vector Ig = ^/||^||. 


3.3. Numerical discretization of the averaging scheme 


Numerically, the helds (14)-(17) can be discretized by replacing the spatial points r with 


a hnite set of points r* which have the coordinates {xi,yi,Zi). On a regular grid each 
point Vi is then the center of a sub-volume V) = AxAyAz which dehnes the averaging 
sub-domain Qp 

Oj = [xi — I Ax, Xi -f \Ax'\ X ... X [zi — ^Az, Zi -1- \ Ax'\ (18) 

and determines the resolution of the continuum representation. As a consequence of 
this coarse graining, it is well possible that inside a sub-domain Dj line segments “cross 
each other” without intersections: although these lines are located on different physical 
glide planes they can not be resolved separately after coarse graining The discretized 


version of the scalar dislocation density (14) is: 


Pi = p ( a ) = ^ 


^ ^ JciNkfi'^k 

k, c(ufc)eOi 


(19) 


In (19) we compute the line length contribution of each dislocation segment by numerical 


quadrature: letting Uk and Wk be the abscissas and weights of the quadrature method of 
choic^in the unit interval [0,1], respectively, each quadrature point contributes a line 
length J{uk)wk to the domain D* which contains c{uk). Conveniently, the quantity J{uk) 
be directly computed for each dislocation segment from (|^ and (§. Analogously, 


can 


f The coarse graining resolution is the reason why short-range interactions have to be handled 
differently as compared to DDD. This, however, will not be discussed here. 

f In this work we have used a uniform distribution of the quadrature points i the interval [0,1]. 
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we define the discretized vector of signed GND densities, and the discretized signed 
cnrvatnre density as 


Qi = Qiri) = Jciuk)wk 


, fc, c{Uk)&fli 


( 20 ) 


<(1 = ‘i(n) = ^ 


E 

. k, c(uk)GQi ^ 


dP 


dc 


X n 


Jcid^k^'^k 


( 21 ) 


Uk 


4. Evolution equations for CDD and DDD 

4-.1. Continuum Dislocation Dynamics 

The Continuum Dislocation Dynamics theory provides conceptual steps to map a 
discrete dislocation system onto a set of continuous, density-like field variables together 
with the respective evolution equations governing their change in time. The original, 
higher-dimensional theory (hdCDD) was defined in a configuration space which added 
an additional dimension to the spatial domain 32,48 . hdCDD is very accurate and 


contains very detailed microstructure information 39,49 but suffers from requiring a 


large number of computational degrees of freedoms. As a computationally favorable 
albeit somewhat less accurate theory a simplified version (CDD) is used, which is 
based on spatial evolution equations for the total dislocation density p, a vector of 


geometrically necessary signed densities g and the curvature density q 35,45,50 . We 


assume that the components of g are oriented such that they represent the orientation of 
screw and edge dislocations: g= [ps 5 ^?e]• The temporal evolution of these field quantities 
in local slip system coordinates - governed by the following set of equations: 

dtp = - V ■ (vg-^) + vq (22) 


are 


dtg = — V X (vpn) 
dtq = -V-(-nQ(^) + A(2).Vn), 

where n denotes the slip plane normal. The by 90^ 
is g^ = [pe,—&], and we assume = —g^q/p (see 
assumption). The evolution equations are mathematically closed by an expression for 
the tensor which here is obtained from a “maximum entropy” approach 


(23) 

(24) 

rotated CND density vector 
for a discussion of this 


50 


50 


AP) = (( [(1 + 4>)C ® C + (1 - 4) V ® 


(26) 


where is the unit vector perpendicular to Ig. As shown by Monavari et ah [50 


$ can be approximated as ~ (p/p)^(l -|- (p/p)^)/2 and is a non-linear interpolation 
between p = 0 (isotropic dislocation arrangement) and g = p (fully polarized dislocation 
arrangement). In section]^ we numerically investigate the validity of this approximation. 


The corresponding non-dimensionless scaling of the CDD equations (22)-(24) can be 


found in Appendix D of 51 . 
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In general, the velocity v in the above equations needs to be specihed in terms 
of a constitutive equation through which elastic interactions enter the model. This is 
a priori not part of CDD: e.g. long-range (or “internal”) interaction stresses have to 


be obtained from the additional solution of a dislocation eigenstrain problem 42,52,53 


while elastic short-range interactions have to be recovered in an alternative way based on 
CDD fields (see e.g. 54 -56 for steps into this direction) due to the limited resolution of 
any continuous description. Within the present work we will not tackle this problem of 
“dynamic closure” and assume for benchmark purposes velocity fields given as stationary 
analytical functions that do not depend on the dislocation stat^ In other words: 
we ignore all dislocation interactions and concentrate exclusively on the “kinematics”, 
that is, how curved and connected lines move and evolve in space and time without 
interactions. 


4-2. Regularization 


Eqs. (19)-(21) can be used to extract coarse grained fields of geometrical data from a 
given configuration of discrete dislocations. However, the CDD evolution equations 
contain spatial derivatives which have to be approximated numerically and which 
thus require field data with a certain level of spatial smoothness. It is therefore 
helpful to replace the Dirac delta-function 5{r) in the convolution integral ([^ with 
a regularization function G{r). To improve the numerical efficiency during computation 
of the convolution it is beneficial to make use of the sifting property of the Dirac function: 
we can always write the regularized function as G = where the symbol * indicates 
convolution over three-dimensional space. Together with the fact that convolution 
and volume averaging commute, it follows that regularized fields can be obtained by 
convolution of (19)-(21) with G. In our numerical implementation, we chose a Gaussian 
standard distribution as the regularization function: 

1 ^ 


G(r) = 




exp 


2s2 


(26) 


where the standard deviation s characterizes the width of the dislocation density 
distribution. For numerical reasons G is approximated by a discrete Gauss function 
G'^ which is non-zero only at a finite number of discrete points r* and is normalized 
such that AxAyAzJ2iG‘^{i) = G(r(x', y', z'))dx' dy' dz' within a small numerical 

tolerance. We then compute the convolution of the coarse grained fields (eqns. (19)- 
(21)) with the discrete Gauss function G^. We note, that in (26) the only free parameter 
is the standard deviation s which subsequently will be chosen as the average dislocation 
spacing s (the mean separation between a dislocation and its nearest neighbors) inside 
the averaging volume (see also the discussion in (^). Note that choosing a much larger 


§ Note that this assumption can be physical e.g. in the case of a high purity semiconductor single 
crystal where the Peierls barriers are high and dislocation densities are low such that the influence of 
dislocation interactions on their motion can be neglected. In general, of course, interactions cannot be 
neglected. 
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value for s would destroy details of the dislocation structure to a large extend while if 
s were very small this would result in a quasi-discrete rather than continuous density 
structure. This is further discussed and demonstrated in the appendix. 

5. Numerical experiments 

One of the key advantages of the CDD theory is the property of kinematic consistency, 
i.e. the ability to represent flowing dislocations through their characteristic geometrical 
properties: a dislocation is a linear, curved line where each of its points moves 
perpendicular to the respective local line tangents. In the following benchmark 
systems we will directly compare the evolution of discrete and continuous dislocation 
microstructure in analytically given velocity helds. We start with 3 situations 
of increasing complexity with regards to boundary conditions: dislocation loop 
distributions in a rectangular domain (i) with open boundaries allowing out-flux of 
dislocations, (ii) with impenetrable surfaces and (hi) with a circular obstacle (i.e. 
an internal boundary). All 3 systems will start from the same initial dislocation 
microstructure which makes it easy to compare the effects caused by different external 
and internal boundaries. The 4th benchmark system is a periodic system with a 
statistically homogeneous random distribution of dislocations and investigates how good 
fluctuation can be represented in a situation without the strong polarizing effect of 
boundaries. 


5.1. System and initial values 


The following examples study a system as sketched in Fig. where a 3-dimensional pillar 
with a square base is discretized into thin sub-volumes along the 2 : axis. We investigate 
one of such sub-volumes with size = Ly = 20006 and height Az = 756, where 
6 = 0.256 nm is the norm of the Burgers vector. Initial values for DDD simulations for 
Systems 1-3 are obtained by randomly distributing Ah = 50 circular, discrete dislocation 
loops with radius R = 756 inside the volume V (note, that we only chose circular 
loops for ease of implementation; loops could be arbitrarily shaped - as long as they 
are closed). Each loop is positioned on a different glide plane and thus loops do not 
intersect by construction; for the hrst three systems their centers are distributed on a 
quadratic domain O.dLa, x 0.5Lj, such that all loops are completely contained inside the 
computational domain. The average density in this volume is po = 4.8 x 10^®m“^. Az 
is the coarse graining height such that the system can be numerically represented as a 
two-dimensional object - a thin lamella in the x — y plane (Fig. [^). Continuous held 
data is obtained by ’smearing-out’ the coarse grained helds through the regularization 


outlined in section 4.2 Initial values for System 4 are obtained by distributing the loop 


centers on the whole quadratic domain x Ly while simultaneously taking care of 
the periodicity. For determining the standard deviation for the Gaussian distribution 
we estimate the mean dislocation spacing from the density pq, which results in a mean 
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dislocation spacing of s = Lxl{2\/2jrNd) ^ L^/Sb. For onr examples this estimate 
of the standard deviation works well for all investigated systems and at all time steps 
(farther details on the choise of s can be fonnd in the appendix). Only when the mean 
dislocation spacing changes drastically - in section [5^ dne to dislocation pile nps at grain 
bonndaries - we adjust this value locally during the post-processing step for obtaining 
reference DDD held data to the actual value of the dislocation spacing. An example 




curvature density q [po/R] 



screw/ GNDs Qs [po] edge GNDs pe [po] 


Figure 2: Initial values for the DDD simulation (upper left) from which initial values 
for the CDD held variables are obtained. 


for DDD and resulting CDD initial structure that is used for all numerical examples is 
shown in Fig. the total dislocation density has a maximum in the center region where 
most of the loops overlap. At the same time nearly equal numbers of loop segments of 
all possible orientations are present in the center region which makes the dislocations 
’statistically stored’. Hence, the GND density |^| tends to zero. The curvature density 
q has the same shape as the total density because initially all loops have the same radius 
and thus q = p/R. Although for each of the hrst 3 studies identical initial values are 
used, we note that the chosen CDD fields are not special or more suitable than any 
other set of statistically equivalent initial values. 

5.2. Numerical methods 

For the numerical solution of the continuum evolution equations a Galerkin finite 
element scheme together with an implicit time integration scheme is used. To accurately 
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represent derivatives of up to second order, Lagrangian shape functions with cubic 
polynomials were chosen. The hnite element mesh consists of 30x30 elements and thus 
has an element size which is in the same range as Az. When high density gradients 
at boundaries occur (System 2 and 3) we additionally rehne elements in those regions. 
To increase numerical stability in particular in regions with high fluxes and where the 
density approaches zero a small viscous term was added to each evolution equation. E.g. 
in ( 22 ) we replace the divergence term V ■ {vq) with V ■ {vg + eVg), where e is chosen 
sufficiently small such that the physical behavior of the equations is not affected while 
numerical oscillations are suppressed. The FEM simulations presented below all take 
less than two minutes computational time on a off-the-shelve computer with a 2.6 GHz 
dual core processor. 

Computation of the evolution of plastic strain can be done for CDD based on the 
Orowan equation {d'jt = pbv) by time integration alongside with the CDD evolution 
equations (22)-(24). To obtain the initial values for 7 from DDD one can compute 7 
from the GND vector and the relation g = geometrical construction of the 

swept area for each unit of GND density and subsequent superposition of these areas - 
followed by division by A 2 ; - gives the plastic strain. Since the plastic strain is a derived 
quantity it suffices to show that main CDD held quantities match those from DDD and 
we do not explicitly show 7 . 

The DDD model used the representation by chordal splines {a = 1); for further 
numerical details on the used Parametric Dislocation Dynamics method please refer 
to 


14,46 


5.3. System 1: distribution of dislocation loops, open boundaries 

For this system a constant velocity of Vg = O.Olnm/ps is prescribed, and time is 
expressed as multiple of the duration that density needs to completely traverse the 
domain, L^/vg. Boundaries are ’open’, i.e. dislocations may leave the volume as if the 
domain were just a sub-domain of a much larger domain. This case of dislocation loops 
which have the same curvature and expand with constant velocity is a case for which the 
CDD theory is supposed to give mathematically exact results, and the system therefore 
can be used as a reference system for testing the quality of the numerical scheme and if 
the chosen standard deviation s is appropriate. Initial values are shown in Fig. Two 
different simulations are run: discrete microstructure is evolved with the DDD model 
and - separately - the continuous microstructure is evolved with the CDD model. To 
compare the two simulations at time t the discrete microstructure was converted into 
CDD held quantities using the same strategy as for obtaining initial values from discrete 
data (in the sketch Fig. which is then compared with E). Fig. shows a snapshot 
in time, where in the contour plots the dashed lines indicate CDD values and the full 
lines are the contours of the converted DDD values. We observe that as dislocation 
loops expand they how away from the center region where the total density p tends to 
zero (compare p in Fig. |^. The diherence in the centers of the loops becomes more 
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System 1 (time t « OAL^/vq) 

1 4.00 

3.20 

2.40 

1.60 

0.80 
Jo.oo 

total density p [po] GND density |^| [po] 




1 2.00 

1.60 

1.20 

0.80 

0.40 

Jo.oo 



screw GNDs Qs [po] 


edge GNDs pe [po] 


4.00 


3.11 


0.44 


- 3.11 


- 4.00 


Figure 3: A snapshot in time (after segments traveled ~ OAL^) for System 1 with 
open boundaries. The color values and solid contour lines show reference DDD data, 
the dashed lines show contours of the CDD data. 


insignihcant as their radii grow and the resulting density distribution takes a loop-like 
shape. In each point dislocations now have all approximately the same orientation, they 
become GNDs which shows in p as can be observed in Fig. The comparison of 
DDD and CDD helds shows excellent agreement. Only minor features of the DDD helds 
are slightly smoothed out which is an artifact of the numerical solution. The volume 
integrals of p and q/2'K give the total line length and the total number of dislocation loops 
and are shown in Fig. For reference we also show the line length for an inhnite size 
system (approximated by a signihcantly larger domain in x and y direction) for which 
the total line length can be obtained analytically as L{t) = NA2 tt{R + vt). Initially, 
both systems exhibit a linear line length increase (loops expand with constant velocity). 
As soon as density leaves the volume through the boundaries the line length for the 
hnite size system decreases; the line length in the inhnite system, however, continues to 
increase linearly. CDD data matches the DDD data perfectly and also coincides for the 
large system with the analytical solution. The right plot of Fig. |^shows the total number 
of dislocation loops. This value stays constant for the inhnite system, which is correct 
since there are neither sources nor annihilation. For the hnite system the surfaces act 
as sinks and the number of loops is reduced. Note that this interplay between density 
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DDD, L,. 


DDD, =2/im, open BCs ▼-T analyt. free loop expansion 




Figure 4: Normalized dislocation line length (with Lq the initial line length) and total 
number of loops for System 1 as function of time. For comparison also the data for an 
unbounded system is shown together with the analytical solution. 


and curvature density will also become relevant when sources and annihilation are to 
be modeled (which we do not consider within the present work). 


5 . 4 . System 2: distribution of circular loops in a finite domain 


Now the outflow of dislocations will be constrained by imposing impenetrable boundary 
conditions, i.e. the dislocation flux pv is enforced to be zero at the boundaries. 
Numerically, this is done by prescribing a velocity held which decays smoothly to zero 
directly at the boundary (in physical terms this boundary layer could be imagined as 
the result of FIB damage in the surface-near regions). To mathematically model the 
boundary layer we use the following sigmoidal function. 


/(O 


l + expK/3] 


(27) 


where the parameter (3 controls the shape of the function and is chosen such that 
the boundary layer is clearly visible and allows to study the dislocation hux in detail 
{(3 = 20). f is the normalized distance from the boundary, e.g. with ^ = 0... 1. The 
velocity held can now be composed from the four contributions 


v{x, y)=vo- f{x/L^) ■ /(I - x/L^) ■ f{y/Ly) ■ /(I - y/Ly). (28) 

It is depicted in Fig. and will be used for CDD as well as for DDD simulations. 
Both simulations are conducted with the initial dislocation structure that was used 
as well for System 1 (Fig. [^. Again, the DDD simulation starts from the discrete 
microstructure, the CDD simulation from the continuous microstructure as obtained 
by D2C conversion. Both simulations evolve independently. Only in a post-processing 
step we again convert the resulting DDD microstructure for reference purposes into 
continuous held quantities as outlined before. A snapshot in time of the respective 
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1 0.010 

0.008 

0.006 

0.004 

0.002 
-lo.ooo 

Figure 5: Prescribed velocity field with uq = 0.01 nm/|us in the center region and n = 0 
at the boundaries of a ’grain’ with impenetrable boundary conditions (System 2). This 
velocity field is used both for the CDD and the DDD simulation 



evolved dislocation field quantities is shown Fig. [^with CDD values shown as dashed 
contour lines and converted DDD field values shown as solid lines (also see the figure 
in the appendix and the supplementary movie 1 which additionally shows the evolution 
of average curvature k = q/p). During the simulation dislocation loops initially expand 

System 2 (time t 
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GND density |^| [po] 
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screw GNDs [po] 
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curvature density q [po/R] 

Figure 6: Snapshot of dislocation structure in a quadratic grain with impenetrable 
boundaries. The solid contour lines show reference data obtained from DDD, the dashed 
lines show CDD data. 
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freely until parts of them reach the boundary layer within which the velocity decays to 
zero. The non-zero velocity gradient rotates line segments such that they eventually are 
parallel to the boundary, i.e. their line tangents are aligned with the contours of the 
velocity held. This causes the initial SSD density to increasingly become geometrically 
necessary when dislocations pile up against the boundaries (compare p and |^| in Fig.[^. 
Additionally, dislocations need to bend with a high curvature radius near the corners 
of the domain. This shows in the high value of the curvature density in these regions. 
For a more detailed view Fig. shows horizontal and diagonal line plots through the 
domain for two different time steps. We observe again an excellent agreement between 
the DDD and CDD held values for all times, and only small deviations can be seen close 
to the boundaries and in particular for the curvature density. 


total density p [po] 




curvature density q [pq/ 



Figure 7: Snapshots in time for System 2: total density (left column) and curvature 
density prohles (right column) for an early point of time (top row) and a later point of 
time (bottom row) which corresponds to Fig. The plotted data is along a horizontal 
line at p = Lyl2 (blue) and a diagonal line ai y = x (red), where the lower left corner 
of the domain is at x = p = 0 (cf. Fig. Q. 


The total line length inside the system serves as a good plausibility test: initially 
we again would expect free loop expansion until hnally all loops are aligned with the 
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Figure 8: Dislocation line length for System 2 as function of time. For comparison also 
the two analytical asymptotes are shown. 


boundary and are (nearly) rectangular in shape. Fig. shows the total line length vs 
time, and it is obvious that at all times DDD and CDD agree very well: both approach 
the horizontal tangent when all lines are aligned with the boundaries. The fact that DDD 
values are slightly below the horizontal tangent is correct because loops are not exactly 
rectangular but still have small rounded corners. For the DDD curve we directly used 
the line length from the discrete data while for the CDD data we hrst had to integrate p. 
The small difference in the hnal value of the line length arises from discretization errors 
due to very steep gradients at the boundaries. In this respect, more realistic systems 
that also consider elastic dislocation interactions will numerically behave even better 
because pile ups always have a hnite size. 

5.5. System 3: circular dislocation loops in a domain with an impenetrable inclusion 

The third system studies the idealization of a plastically deforming matrix into which an 
elastic, circular inclusion is embedded. Dislocations can move freely in the matrix while 
the inclusion acts as an impenetrable obstacle for dislocation motion. Annihilation of 
lines is not considered, hence, the number of dislocation loops stays constant; for real 
systems this is a strong restriction, but again, the philosophy is to decouple interactions 
and reactions from the purely kinematic effects. The computational model is - as in 
System 1 - a quadratic domain = Ly) with open boundaries and contains the 
inclusion in the center. The inclusion itself is represented by an internal boundary 
which we model as a circular sub-domain (radius R = L^/IO) with zero velocity in the 
center region and a smooth transition from n = 0 to a constant velocity uq away from 
the inclusion (cf. Fig. |^a)). As a result of the velocity held dislocations are hrst slowed 
down on their way towards the inclusion and eventually freeze in their motion so that 
they cannot enter the inclusion. Both simulation methods use the same initial values 
as before. This does place some dislocations inside the inclusion which is unphysical. 
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velocity v [nm/ [xs] 



Line length vs time compared to System 1 


Figure 9: Velocity field with the “inclusion” in the center and evolution of line length 
(normalized with initial total line lenght Lq). 


It turned out, however, that the influence on the resulting microstructure is only small. 
Therefore, we decided to consistently use the same initial values as for all other systems 
as well. Evolving the DDD and CDD system gives typical dislocation patterns as shown 
in Fig. 10, additional plots including error plots can be found in the appendix. 

An auspicious feature is the formation of pile-ups of bent lines around the inclusion: 
when expanding dislocation loops touch the inclusion they get a ’dent’ (see the colored 
line in Fig. [Io|(a)) and adjust to the shape of the inclusion. At the same time, line 
segments further away start ’spiraling’ around the inclusion and deposit an increasing 
number of segments around the inclusion (also see the hgure in the appendix and the 
supplementary movie). This behavior has a number of consequences and can be observed 
in the CDD representation as well: 


• Circular pile-ups of dislocation density can be seen in the upper middle and right 
plots of Fig. [T^ : p has a high value in the vicinity of the inclusion because density 
from piled-up dislocations together with pre-existing line segments in the inside of 
the precipitate add up. In conjunction with the GND density |^| we infer that 
dislocations in the inside are mainly statistically stored while pile-ups around the 
inclusion are - literally - geometrically necessary. 


Circular CNDs around the inclusion form an inverted loop, cf. Fig. 10 bottom right 
and middle plot. Taking a look at the sign of e.g. the edge GNDs Pe around the 
inclusion it is found that their orientation was rotated by 180° (also compare the 
initial values in Fig. |^. For instance the upper half is now negative and therefore 
would move downwards but is hindered by the inclusion. 


Bent dislocations around the inclusion have negative curvature - the above 
statement about an inverted loop already suggests that the curvature of dislocations 
around the inclusion should be negative which in fact can be observed in the two 
plots on the right of Fig. 11 where the average curvature k = q/p was computed and 
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System 3 (time t « Q.iL^/'vf) 



DDD total density p [po] GND density |^| [po] 



curvature density q \po/R] 



Figure 10: Evolved dislocation loops in a domain with an precipitate, i.e. a region in 
the center with zero velocity. The color values and solid lines show reference DDD data, 
the dashed contour lines are CDD data. 


plotted along two lines passing through the inclusion. The inclusion has a radius of 
Lj-ZlO which is equivalent to a curvature of fcjnc = 0.020nm~^. Due to the smooth 
interface the velocity is zero only at about L^/IQ which is equivalent to a curvature 
of fc min = 0.031 nm“^. In the bottom right plot of Fig. 11 we find the negative 


curvature peak values at |A:cdd| ~ 0.040 nm ^ in good agreement with both the 
geometry value /cmin and the DDD curvature value of IfcoDol ~ 0.028 nm 


-1 


Deviations between DDD and CDD are larger than for the previous two benchmark 
models. They are particularly large in those regions which exhibit large gradients in 
the velocity and in the CDD held values, both of which cause increasing discretization 
errors. What is the average response of this system? Comparing the temporal evolution 
of the DDD line length in Fig. IT to System 1 shows that the inclusion creates additional 
density which is clearly visible from approximately t = t).lbLy;/vQ on and is caused by 
lines bending around the inclusion. CDD is lagging behind and is initially even slightly 
below the line length of System 1 (due to the lower average velocity as compared to 
System 1). After t > t).2'bLx/vQ the CDD model is then also showing an additional 
increase in line length. 

The difference between CDD and DDD is the result of larger errors close to the 
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total density p [po] curvature k = q/p [nm 




Figure 11: Snapshots in time for System 3: total density (left column) and curvature 
profiles (right column) for an early point of time (top row) and a later point of time 
(bottom row) which corresponds to Fig. 10 The plotted data is along a horizontal line 
at p = Ty/2 (blue) and a diagonal line a.t y = x (red), where the lower left corner of the 
domain is at x = p = 0 (cf. Fig. Ill) 


inclusion where CDD almost always underestimates the true density values. Increasing 
the resolution of the spatial discretization of the hnite element scheme does not yield any 
appreciable alleviation. The reason is the very complex deformation state of dislocations 
that brings the CDD theory to its limits because near the inclusion the line curvature 
can be different for different line orientations in the same averaging volume - a detail 
that cannot be represented with this simplihed variant of CDD (this CDD theory only 
can represent one average curvature value for each point/averaging volume). Although 
a number of details of this system have been predicted properly - e.g. the line curvature 
around the inclusion which is important when it comes to hardening effects in terms 
of line tension - this particular CDD formulation only roughly predicts the line length 
production close to the inclusion. More elaborate evolution equations which contain 
more information about the orientation distribution and curvature of dislocations are 
required as e.g. those derived in 45,50 . The D2C strategy for more detailed CDD 
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formulations, however, will still remain the same. 

5.6. System f: random dislocation distribution in a periodic domain 

The previous benchmark systems have in common that they start with random 
initial values which during time evolution become strongly polarized through the 
geometrical constraints imposed by boundary conditions (e.g. dislocations adjust their 
line orientation and curvature to the shape of external or internal boundaries and become 
geometrically necessary). System 4 will now investigate how good CDD performs in a 
bulk-like situation where the total density consists to a large extend of statistically stored 
dislocation density superimposed with (smaller) GND density fluctuations. Periodic 
boundary conditions are used to eliminate geometrical constraints. Initial values consist 
of a statistically homogeneous random distribution of 200 dislocation loops of the same 
radius, and their centers are distributed across the whole volume L^. x Ly x Az. This 
results in the same average density po as in the center region of Systems 1-3; all other 
parameters stay the same. 


System 4 (periodic boundary conditions) 



DDD (p - (p)) /var(p) (|e| - (|^|)) /var(|e|) 


Figure 12: Initial values and evolved state of dislocation microstructure in a periodic 
domain. Contour plots show the normalized fluctuations of total density (middle 
column) and of the GND density (right column). The color values and solid lines 
again denote reference DDD data, the dashed contour lines are GDD data. 
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DDD: - — (p) 

CDD: ---- ---- pSSD .... ^GND .... pfluct .... 



Figure 13: Densities along a horizontal line through the center of the computational 
domain for Study 4: initial state (left) and evolved state at t = OAL^/v (right). 


The velocity and loops’ curvature radius are constant everywhere. Since this study 
analyzes the quality of the representation of flow and evolution of density fluctuations, 
it is useful to decompose the total density into fluctuating and volume-averaged 
contributions, p{x,y) = pauct{x,y) -|- (p). For plotting the held data we additionally 
normalize the huctuations with the variance of the respective held, e.g. in Fig. [T^ we 


plot for the total density {p— (p))/var(p). Fig. 13 shows line plots of initial and evolved 
DDD and CDD data, additional error plots can be found in the appendix. Clearly, all 
details of the discrete microstructure are accurately reproduced. Fluctuations evolve in 
a non-trivial manner and we attribute negligible deviations to numerical errors. The 
average behavior matches exactly the DDD values. The overall very good agreement is 
not surprising after seeing that already System 1, which initially had a similar random 
structure in the center of the domain, gave very good agreement between DDD and 
CDD. It is, however, important to show that this good agreement is sustained during 
time evolution and that even the increase in the range of the huctuations are properly 
represented. An application based on this type of system has been studied in the 
context of dislocation patterning 41 where in particular the huctuations together with 
the correct line length increase and a Taylor-type how stress were identihed as key 
ingredients for pattern formation. 


6. Determining the kinematic closure parameter $ from DDD 

The numerical studies in section [5] demonstrated that DDD simulations can serve as 
an elegant means for directly benchmarking continuum dislocation microstructure and 
its evolution. Thus, it seems natural to directly use DDD simulations to verify or 
even obtain expressions that were used to constitutively close the set of CDD evolution 
equations. Within the CDD equations occurred in particular one quantity which had to 
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be assumed: the interpolation function $, which was taken in the ’maximum entropy 
approximation’ approach as = (|^|/p)^(l + (|^|/p)^)/ 2 . $ was in an earlier approach 
assumed to interpolate linearly between SSD and GND density. We will now analyze 
DDD simulations in order to decide which assumption is under realistic simulation 
conditions the most suitable one. Recall the definition of the evolution equation of the 
curvature density and 

_ _ div(-nQ(^) + A(2) . Vn), (29) 

AP) =^[(1 + 4)(,0(,+ (1-4)(^0(P]. (30) 

In A*-^) the approximation occurs in a vector product together with a spatial 
velocity gradient. Thus, dtq^^ is sensitive w.r.t. variations in only in regions where 
the velocity field is non-constant, as e.g. in the boundary layer close to the impenetrable 
boundary of System 2 from section 5A Continuous CDD field variables can be obtained 
from a DDD simulation as outlined before, and the time derivative dtq can be numerically 
approximated by an explicit Euler scheme in time. 


dtq^^ = 


q^^{ti+i) - q^^jtj) 

At 


(31) 


At the same time also the evolution equation for q (eqn. (29)) can be evaluated based 
on field data obtained from DDD and Vn are both known, and I, l'^ are functions 
of g and p)[j|] The unknown value of can now be obtained for each point (where 
Vv 7 ^ 0) from solution of the inverse problem : 

Find <h(x, y) — for all (x, y) with Vv{x, y) ^ 0 — such that 

\\dtq^^{x,y) - dtq^^{^{x,y);x,y)\\ MIN. 


Note that no CDD simulation needs to be done, only DDD data is used. For each 
point {x, y) we record the optimum value of $ together with the local total density 
p{x,y) and CND density g{x,y). Plotting these data gives the distribution shown in 
Fig. There, the point of time was chosen such that already a significant amount 
of density reached the boundary layer but did not yet reach a stationary state. For an 
earlier point of time dislocations would not yet have reached the region where Vn 7 ^ 0 
and no data could be obtained. At a later point of time most dislocations form pile- 
ups of geometrically necessary configurations near the boundaries which results in a 
strong accumulation of data points at the top right corner of the diagram. It is obvious 
that the optimization strategy yielded points with a relatively large scatter. However, 
it can be concluded that most data points are located below the straight dashed line 
which is the linear interpolation between CNDs and SSDs. The polynomial maximum 
entropy approximation fits much better. Binning the data e.g. in 20 g/p intervals gives 
averages that for g/p = 0 . 8 .. 1.0 coincide nicely with the polynomial approximation. 


II We note, that the space dependent tensors and also can be extracted from DDD data. 
E.g. the TP component of can be obtained by integrating the discrete density for each separate 
line segment against the cos^ ipc of the discrete line orientation ipcj followed by regularization based on 
the Gaussian convolution. 
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Nonetheless, the data is not sufficient for any statistical analysis in all other regions 
and could not be much improved through additional simulations. These results might 
indicate that the approximation of dtq is not sufficiently sensitive w.r.t. to the only 
parameter q/ p. Currently, higher order closure assumptions are under investigations. 
Their results might help to understand the difficulties in identifying a unique functional 
$ for closure. 



Figure 14: <h(p/p) for the linear approximation and the ’maximum entropy 

approximation’. Data points show best hts of the numerically computed values for 


7. Conclusion 

We developed a systematic method for averaging geometrical properties of discrete 
dislocation lines which allows for a direct comparison of discrete and continuum 
descriptions of evolving dislocation microstructures. Particular care was taken to 
formulate the numerical approximation such that discrete and continuous formulations 
are at all times consistent with each other (e.g. both the DDD and the CDD 
descriptions were based on the same spline representation). Using DDD data as reference 
the comparison of simulation results for 4 different systems revealed a surprising 
accuracy of the CDD theory and its numerical implementation: CDD is able to evolve 
very complex dislocation microstructure and simultaneously respect various physical 
boundary conditions in almost exactly the same level of detail as DDD simulations 
do. Even in situations where the simplifying assumptions of the used CDD formulation 
were clearly violated (System 3) we found that results were still reasonably. Within this 
work we only considered the kinematic closure and assumed a (hxed and analytical) 
velocity function. An important step for future work which still requires a large amount 
of fundamental work, however, is the “dynamical closure”, i.e. the question what the 
average velocity of interacting continuum distributions of dislocations under stress is. 
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Considering that in 3D DDD simulations the computational cost scales 
approximately with {N being the number of nodes/segments) and considering that 
for a continuum theory the computational cost does not increase with the number 
of dislocations at all (not even when dislocation interactions are considered) we are 
very optimistic that a continuum description of dislocation plasticity might very soon 
complement DDD simulations in particular when it comes to high densities in large 
volumes and/or high accumulated plastic strains. Furthermore, the D2C conversion 
might even be directly beneficial for analyzing DDD data: the continuum representation 
is well suitable for ensemble averaging and could be a novel way of directly comparing 
discrete microstructures. 
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Appendix. On the resolution of continuum fields 


In section continuous fields were obtained that had to be artificially smoothed out 
for numerical reasons. This was done by convolving the coarse grained fields (eqns. 
(19)-(21)) with a discrete Gauss function Ghoosing the standard deviation s in (26) 
is a choice which has to be made for every system and e.g. with regards to numerical 
aspects as well as based on the desired degree of detail which the density fields should be 
able to represent (e.g. in terms of density fluctuations). In the following, we summarize 
some key aspects that one might want to consider for deciding on a suitable value for s. 


• The obvious criterion for a lower bound for s is related to the resolution of the 
numerical scheme: choosing s much smaller than roughly the size of a finite element 
(which additionally depends on the used shape functions) results in discretization 
errors, i.e. the fluctuations of the density field cannot be numerically represented 
properly anymore. In principle the discretization error can be measured. For 
practical purposes one could e.g. compare the numerical representation of a 
’continuum loop’ with the analytical solution. Note that choosing a value of the 
standard deviation in the range of a Burgers vector results in a computationally 
very expansive continuum model of nearly discrete objects (Fig. |Al|left). 


A numerically stable solution of the transport equations of density fluxes 
additionally requires a certain degree of smoothness in order to avoid strong, 
undesired numerical oscillations. Studying the aspect of numerical stability is non¬ 
trivial (see e.g. 57 for a stability analysis for a continuum dislocation dynamics 


model based on the Kroner-Nye tensor) and was not attempted for the present 
GDD equations so far. We approach this problem in a pragmatic way and choose 
numerical parameters as well as the lower bound for the standard deviation s 
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such that a single continuum dislocation loop can be expanded accurately within a 


small error tolerance during the simulated time 48 


Deciding on an upper bound for the standard deviation s is problem specihc: 
CDD works properly for all values of s exceeding the above introduced lower 
bounds. However, the larger the chosen s the more details of the heterogeneous 
microstructure will be destroyed. In particular, choosing s in the range of magnitude 
of the system dimension causes all gradients of the density helds to vanish. As a 
consequence, the partial differential equations for transport of densities deteriorate 
to ordinary differential equations, which are no longer able to describe fluxes on a 
scale below the system size (Fig. Al right). 


In this work, the emphasis is on heterogeneous dislocation microstructure. Hence, 
we chose, guided by the mean dislocation spacing of s ~ L^/35, a value of s which is 
small enough that fluctuations are not getting smeared out, but which at the same time 
is large enough so that discrete dislocations can no longer be differentiated. 


s = Lj-ZdOO s = L^/AQ s = L^/A 

1 4,95 

4,12 

3,30 

2,47 

- 1.65 

- 0.82 
-J-0.00 

Figure Al: Initial density p for different values of the standard deviation s obtained 
from the DDD data of System 4. The left plot is a nearly discrete microstructure, 
while in the right plot an almost homogeneous structure without any gradients resulted. 
CDD can in principle evolve all initial values accurately, the only limitation is our used 
numerical scheme, which does not work for values of below s ~ 




Appendix. Additional time evolution and error plots for System 2, 3 and 4 


In Fig. B2 - B4 the time evolution of discrete and continuous microstructure together 
with the relative error — p ^^\/for System 2 (grain with impenetrable 

boundaries). System 3 (’precipitate’) and System 4 (periodic BCs) are shown. In order 
to reduce the unreasonable diverging behavior of the relative error in regions where 
—)■ 0 we chose a cut off value of 0.5 (i.e. less than half the size of the smallest 
interval of the color scale of the density contour plot for System 2 and 3) below which 
the error is not computed (shown as white regions). In the density plots the solid lines 
are - as before - the contours of the converted DDD data while the dashed lines show 
the CDD data. 
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System 2 in Fig. |B2| has in most regions and for all time steps relative errors of less 
than 5%. Only at later time steps when dislocations pile up against the boundaries 
and density exhibits steep gradients deviations become larger in some places. The CDD 
theory does represent the evolution of this system properly. System 3 in Fig. |B3| shows 
clearly that the geometrical constraint of the inclusion creates a complex situation which 
only partially can be represented by CDD and causes large errors. Despite the fact that 
a number of features are correctly represented (in particular the curvature close to the 
inclusion, see main text), a more precise evolution of this system would require a more 
rehned CDD theory. System 4 in Fig. BA demonstrates that CDD consistently and 
accurately predicts the increase of average density as well as the transport of density 
fluctuations - both of which are automatically contained within the set of evolution 
equations. The relative error is in most places even smaller than 4%. 



Figure B2: System 2 - DDD data, density and relative error at time steps f = 0, 0.2 
and O.ALx/vo (from top to bottom). The white regions in the error plots indicate that 
the relative error was not computed due to (nearly) zero density. 
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Figure B3: System 3 - DDD data, density and relative error at time steps t = 0, 0.135 
and 0.27Lx/vo (from top to bottom). The white regions in the error plots indicate that 
the relative error was not computed due to (nearly) zero density. 
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Figure B4: System 4 - DDD data, density and relative error at time steps t = 0, 0.25 
and 0.5L^/vo (from top to bottom). Note the changing color range for the total density. 
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